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Abstract 

We discuss the roles of continuum linear elasticity and atomistic cal- 
culations in determining the formation volume and the strain energy of 
formation of a point defect in a crystal. Our considerations bear special 
relevance to defect formation under stress. The elasticity treatment is 
based on the Green's function solution for a center of contraction or ex- 
pansion in an anisotropic solid. It makes possible the precise definition 
of a formation volume tensor and leads to an extension of Eshelby's re- 
sult for the work done by an external stress during the transformation 
of a continuum inclusion {Proc. Roy. Soc. Land. Ser. A, 241 (1226), 
376, 1957). Parameters necessary for a complete continuum calculation 
of elastic fields around a point defect are obtained by comparing with an 
atomistic solution in the far field. However, an elasticity result makes it 
possible to test the validity of the formation volume that is obtained via 
atomistic calculations under various boundary conditions. It also yields 
the correction term for formation volume calculated under these boundary 
conditions. Using two types of boundary conditions commonly employed 
in atomistic calculations, a comparison is also made of the strain energies 
of formation predicted by continuum elasticity and atomistic calculations. 
The limitations of the continuum linear elastic treatment are revealed by 
comparing with atomistic calculations of the formation volume and strain 
energies of small crystals enclosing point defects. 

* Assistant professor, Dept. of Mechanical Engineering, and Program in Applied Physics. 
Corresponding author. krishnaSumich . edu 

t Assistant professor, Dept. of Materials Science and Engineering, and Program in Applied 
Physics. 

■^Formerly, graduate research assistant. Dept. of Materials Science and Engineering. Now 
a post-doctoral fellow at Institute of Materials Research and Engineering, Singapore. 
^Graduate research assistant. Dept. of Materials Science and Engineering. 
^Graduate research assistant. Dept. of Mechanical Engineering. 



1 



1 Introduction: Free energy of point defect for- 
mation under external stress 



When a point defect forms in a crystal there are two main contributions to the 
Gibbs free energy of formation, = ("^^ — TS^^) — y^™. The first comes 
from the internal energy of formation, and the usually negligible entropic 
term —TS^^, where T is the temperature and 5^^ is the change in entropy. The 
second contribution arises due to the work done by the external stress, and can 
be formally written as —W°^ = —<t°^: V^, where cr'^^ is the external stress 
and is the "formation volume tensor" . The latter quantity is commonly- 
used in the materials physics literature (see Aziz, 1997; Zhao et al., 1999a, b; 
Daw et al., 2001) and can be formally defined as = —d^S'^jdcF^^. A related 
quantity, the migration volume tensor, V" = —d'^^jdcr'^^^ is associated with 
the free energy difference the ground state and the transition state that the 
system traverses when the defect hops between lattice sites. These quantities 
are of fundamental interest in materials physics. They influence the kinetics 
of transport through the diffusivity, D — I?oexp[— (^^^ + ^^™)/fcBr], where fc^ 
is the Boltzmann constant.^ The thermodynamics of transport is influenced 
through the chemical potential, [i = [1q — cr'^^: (V^ + V"), where /io is the 
potential without stress effects. 

It is natural to attempt a continuum treatment of point defects since their 
interaction with stress — a continuum quantity — is of central interest in this pa- 
per. Point defects are modelled as centers of contraction or expansion within 
continuum linear elasticity. While the idea of a formation volume tensor, intro- 
duced above in association with a point defect, is not standard in continuum 
mechanics, we place this notion on firm footing in Section 3. Additionally, the 
introduction of a formation volume that is energy-conjugate to the stress does 
have an analogue in linear elasticity: In a series of seminal papers on elastic 
inclusions Eshelby considered "cutting and welding" operations to describe the 
transformation of an inclusion within a solid that is itself under external stress 
(see, for instance, Eshelby, 1957, 1961). He was able to show that the work 
done is / cr'^^ : e^<W ^ where is the transformation strain due to relaxation, 
and the integral is over the inclusion. We also demonstrate in Section 3 that 
an extension of Eshelby's result to point defects can be obtained in a rigorous 
manner, leading to the form ~W°^ = —(y°^: introduced above. 

The other goal of this communication is to discuss practical aspects sur- 
rounding the evaluation of and . On modelling a point defect as a center 
of contraction or expansion one can obtain a deformation field that clearly must 
bear some relation to and 'i/^ . Whatever these relations, and cannot 
be calculated unless the strength of the center of contraction or expansion, rep- 
resented by a force dipole tensor, is specified in the elasticity problem. Indeed, 
in linear elasticity, all elastic fields scale with the dipole tensor. This quantity 

^Formation plays a role in self-diffusion of point defects, and in inter-diffusion, since both 
phenomena require the formation of a point defect that subsequently migrates. In doing so it 
may enable the migration of a substitutional atom (inter-diffusion). 
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must be obtained from some form of atomistic calculation, as is done in Sec- 
tion 3. However, the linearity of the theory enables us to show in Sections 3 
and 4 that even without quantifying the dipole tensor, rather simple elasticity 
calculations reveal the influence of boundary conditions that are used in the 
atomistic simulations. A straightforward scaling with the correct dipole tensor 
then leads to the actual fields in elasticity. These calculations are carried out for 
the isotropic vacancy in silicon in this preliminary communication. The results 
are summarized and placed in the proper context in Section 5. 

2 The continuum elastic model of a point defect 

In continuum elasticity a point defect is modelled as an arrangement of point 
forces in a center of contraction or expansion. An associated dipole tensor is 
defined as^ 

3 

D = ^Fi®2dj. (1) 

1=1 

In the dipole illustrated in Figure 1 the forces are in equilibrium. Moments are 
also in equilibrium provided that D is symmetric, which will be assumed here. 
In this case D can be written with respect to an orthonormal set of basis vectors 
{61,62,63} that coincide with its own eigen vectors. In terms of this basis, D 
is therefore diagonal. Furthermore, if the magnitudes of the point forces are 
equal: |i^i|, |i^2|, l^^al ~ F, and the magnitudes of the position vectors are also 
equal: |di|, |d2|, jdsl = d, then D is isotropic and can be written as 

D = Dl, where D = ±2Fd, (2) 

and 1 is the second-order isotropic tensor. The positive sign in (2) holds for 
centers of contraction and the negative sign is for centers of expansion as is easily 
verified. The isotropy of D, and D > have been assumed for all numerical 
calculations in the paper. These choices correspond to an isotropic vacancy — 
one in which the magnitudes of displacements of the nearest-neighbor atoms 
are equal, or a substitutional defect with smaller atomic volume than the host 
atoms. If Z? < 0, it is an isotropic interstitial or an isotropic substitutional 
defect^ with larger atomic volume than the host atoms. 

■^Both direct and index notations are used in this paper for brevity and clarity as deemed 
necessary. 

^Anisotropic interstitials, substitutional complexes and even vacancies are possible. In par- 
ticular, one of the configurations yielded by quantum mechanical calculations of the relaxation 
around a vacancy in silicon is a tetragonal distortion of the four nearest-neighbors. In the 
equilibrium position after vacancy formation, all four atoms are at the same distance from the 
vacancy. However, the distances between the four atoms are not the same: The atoms form 
two pairs such that the distance separating the atoms in a pair is smaller than the distance 
separating each atom in one pair from the two in the other pair. This is the Jahn- Teller 
distortion, and in some instances has been calculated to have a lower formation energy than 
the isotropic distortion (Puska et al., 1998; Antonelli et al., 1998). 
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Figure 1: An arrangement of point forces leading to an isotropic dipole — 
specifically an isotropic vacancy or substitutional defect that is smaller than 
the host atoms. 

Using (1), the force distribution due to the dipole is then written as 

3 

fix) = -2Y,F,(E>d,V^,S{x~x') 

i=l 

or fix) = -DV^,6ix-x') 

where S{x — x') is the three-dimensional Dirac-delta function that satisfies 
J hix)5ix — x')dV = hix') for a sufficiently smooth (C°°) test function, hix). 
Since Vj;'(5(a; — x') = —'VxSix — x'), we finally have 

fix)=DVJix-x'). (3) 

The force distribution can be combined with the infinite space Green's func- 
tion for anisotropic linear elasticity to obtain the displacement fields around the 
defect. The Green's function for elasticity, G(a; — x'), is a second-order tensor. 
For an anisotropic, linear medium it satisfies 

Cyfci + S^mSix - x') = 0, (4) 

OXjOXi 

where Ciju is the fourth-order anisotropic elasticity tensor, Gkm is the displace- 
ment at X along due to a unit point force acting at x' along Bm, and 5im 
is the Kronecker delta symbol. Barnett (1972) has applied the Fourier trans- 
form to derive analytic formulae for Gir, its first and second spatial derivatives. 
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These formulae are summarized here: 



G^r = ^^.^ / MU^m^% (5) 



j {-TsMt + ZsH,r) d*, (6) 



dxs Att'^Ix — 



dx sdxjYi 



Letting k denote the wave vector in Fourier space. 2; fc/|fc| is a unit vector 
in Fourier space; M*M = 1, where Mir{z) ^ CijrsZjZg; T = {x — x')/\x — x'\ 
is a unit vector in real space; "if is the polar angle in the plane z - T = 0. Using 
these tensors and vectors H and A are defined as 

A„ = Cjp„^ [{zpT^ + z^Tp) {H,,M*^ + M*^H^r) - 2M*jM*^TpT^] . 



2.1 The displacement and strain fields of a point defect 

The displacement field of a point defect can be written using (3) and the Green's 
function: 

u°"{x) = J G{x - x')DV^,5{x' - x")AV' 



where the superscript on the left hand-side signifies the infinite space solution. 
Observe that the variable of integration is x' . Reverting to indicial notation for 
clarity, and using a standard result on derivatives of distributions, this gives. 



uT{x) 



_d_ 



yGij {x — x') 



D,k5(x' - x")dV'. 



Using dGij {x ~ x')/ dx'^. = —dGij {x — x')/dxk , and the definition of the Dirac- 
delta function we have 



u-(x) 9G,A---") ^^ 
u, [X) - U,k. 

For an isotropic dipole, Djk = DSjk, this simplifies to 



urix)^D—G,,{x-x") 
dxj 



(8) 



(9) 
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The strain field is written as 



1 / a^G,;, d^G, 



D,k, (10) 



2 \dxkdxi ' dxkdxi 
which, for an isotropic dipole, simplifies to 

.oo^j^l 

2 \dxjdxi dxjdxj ■ ^ ' 

Evaluation of (8) and (10) involves the direct application of (6) and (7) respec- 
tively. 



3 The formation volume tensor 
3.1 The influence of boundaries 

The formation volume tensor of a point defect can now be established using 
the center of contraction or expansion model. The results of Section 2.1 for 
infinite crystals are first extended to finite crystals in order to compare with 
the atomistic calculations to be described in Section 3.3. The development in 
this subsection was suggested to the authors by Barnett (2004) for the scalar 
formation volume, tr[V^^]. The extension to the full tensor has not appeared 
before to the knowledge of the authors. 

The formation volume tensor, V^i, is the sum of tensorial volume changes 
due to relaxation of the crystal around the center of contraction or expansion, 
V^i^ and the addition of an atomic volume, ^Vl5ki. The second term arises since 
the displaced atom's volume is added to the crystal. Note the assumption of 
isotropy associated with this step. We obtain an expression for V^^ below that, 
in addition to having other implications, demonstrates that V^^, and therefore 
V^i, are well-defined quantities. 

For a defect in a finite crystal, i?crys, we define 

VL = I eudV + ^nSki, (12) 

-Bcrys 



where Eki is the strain field in the finite crystal due to the center of contraction 
or expansion with the appropriate value of the dipole. Using the stress-strain 
relation, this is 

V^i = SkHj J <7,jdV + ^nSki (13) 

-Bcrys 

where SkUj is the constant compliance tensor satisfying SkUj Cijmn 
fourth-order symmetric identity tensor. Observe that 



d{Xi(J„ij) _ dcTrnj 

Xj 



dV. (14) 



-Bcry 
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Using (3) the equilibrium equation for the center of contraction or expansion is 
written as 



^+dJ-%^ = 0. (15) 
Combining (14) and (15) and using the symmetry of a we have, 



~r XiDjm p. 

KJ X irfl. f-' X fi 



dV. 



Bays 



Using the standard result for spatial derivative of the Dirac-delta function, the 
symmetry of cr, and the Divergence Theorem, 

Bcnys Scuys 

where S'crys is the surface of the crystal. Substituting this equation in (13) we 
have 

V^l = SklijDji + Sklij J Xia-jraUmdA + -^Skl- (16) 

5'ci:vB 



This result has three important, related, consequences: (i) The relaxation 
and formation volume tensors for a crystal with traction-free boundaries depend 
only on the strength of the center of contraction or expansion, represented here 
by the dipole Dji, elasticity of the crystal, and the atomic volume (for the case 
of formation volume). In particular, neither the shape and size of the crystal 
nor the location of the center of contraction or expansion inside the crystal 
influence the outcome, (ii) The difference in V^i between the case with traction- 
free boundary and one with arbitrary boundary conditions can be evaluated. For 
pure traction boundary conditions this is trivial. With displacement boundary 
conditions the elasticity problem must be solved first, possibly by numerical 
methods, (iii) Formation and relaxation volume tensors, unorthodox notions in 
continuum mechanics, can be defined in a precise fashion. 



3.2 A thermodynamic basis for 

The center of contraction or expansion model for a point defect thus leads to a 
relation for that can be evaluated if the stress field of the crystal is known. 
However, the thermodynamic basis of this quantity has not been clarified beyond 
the formal relation = —d'^'^ j dcr'^^ . In this subsection it is demonstrated that 
a result of Eshelby (1957), established for continuum inclusions, also holds in 
the center of contraction or expansion limit and provides such a thermodynamic 
basis. 

We begin with the continuum result for the work of interaction between 
the traction, cr'^'^n, applied at the boundary of a crystal, and the deformation 
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induced by a transforming inclusion in the crystal. As alluded to in Section 1, 
Eshelby showed that this work of interaction has the form 

= y" o-"'': e'dV" (17) 

where i?inc is the region occupied by the inclusion before it undergoes the stress- 
free relaxation strain e'' (Eshelby, 1957). This is the strain in the inclusion if 
its boundary. S-mc, is traction-free. If the inclusion is small enough that cr^^ is 
constant over i?inc, then following (12) and defining Vl^^^ = /g e^'dV allows 
us to rewrite (17) as 

^cx^^cx. yr^^_ (18) 

A seemingly obvious step is to link Vl^^ and V for the point defect. Then 
Eshelby's result would apply to the work of relaxation around a point defect, and 
the work of formation would be a'^^: (V + -jf^l). Introducing this work term 
into the Gibbs free energy would lead to = —d'S'^jdcr'^^. Such an approach 
has been adopted by Aziz et al. (1991). In our notation their expression is = 
^ex . . Vol(-Binc), whcrc is the uniform transformation strain undergone by 
a "sub-system", -Bine- However, there arc technical difficulties in making such a 
direct association: (i) The boundary between the continuum inclusion and the 
solid in which it is embedded is well-defined, but it is not possible to precisely 
demarcate a lattice point defect in this manner. For example, Figure 2 shows 
a vacancy in a silicon crystal where the (10 0) direction is perpendicular to the 
plane of the figure (this is a visualization of one of our atomistic calculations). 
It is evident that no precise interface exists between the defect and the perfect 
lattice. Indeed, in the above expression used by Aziz and co-workers the extent 
of the "sub-system" is not made precise beyond a statement that it contains the 
atoms "involved in the fluctuation" of the transformation. The transformation 
strain is also not made precise, (ii) The arguments that Eshelby used in arriving 
at (17) require that the continuum assumption holds even in the neighborhood 
of the inclusion. The discreteness of the lattice clearly negates this assumption 
in the neighborhood of a point defect. 

The center of contraction or expansion model of the point defect in contin- 
uum linear elasticity enables a circumvention of these difficulties: Consider a 
finite crystal with traction-free boundaries that is also stress-free in its interior. 
As before, i?crys is the domain of the crystal and its boundary is S'crys- Let an 
external traction, cr^^n be applied at S'crys- Also let e'^'^ and u'^^ be the related 
strain and displacement fields in i?crys- Now let a point defect form in the in- 
terior of -Bcrys with corresponding dipole tensor D. Let the stress, strain and 
displacement fields arising from this defect be cr^, and viF in the case that 
5crys is traction-free. Also recall that according to the linear theory of elasticity 
the fields actually obtained now are <t™ -|- cr^, -I- and -|- 'uF . The 
quantity of interest is the work done by the external stress on the crystal in the 
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Figure 2: A vacancy in silicon viewed along the (1 0) direction. The open circles 
mark positions of atoms in the perfect crystal; i.e., before vacancy formation. 
Displacements have been scaled 10 x for clarity. 



process of defect formation: 

Scrys 

Using the Divergence Theorem and integration by parts this is 



5r-= J +4a-)dy, 



where the symmetry of cr^^ has been used to restrict the displacement gradient 
to its symmetric part, the strain. Since (T°^ is divergence-free in i?crys the first 
term in the above integrand vanishes. The second term is rewritten using the 
stress-strain relations to give 



Invoking the symmetry of cr^ and using integration by parts, 

Scrys -^crys 

Since is obtained as the solution to the traction-free boundary case, the 
first term in this integrand vanishes. The second term is rewritten using afj j + 
Dij6j{x — x') = 0, for a defect aX x = x' , to give 



Scrys 



/ <"D,~5{x - x')AV. 
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Using the standard result for derivative of the Dirac-delta function this gives 
For a symmetric dipole tensor (cf Section 2) 

From the stress-strain relations in the form ef^ = ^ijkio'ki and the major and 
minor symmetries of Sijki we have 

Finally, recalling (16) and noting that ajmnm in the boundary integral of that 
equation is now replaced by aj^rim, = we have 

y^'"" '^iMi- (19) 

On including the contribution of the displaced atom as in Section 3.1, this 
result can be extended to the work done in defect formation: 

y^'"" = '^tf (v^i + l^S,i^ = '^uV^i ■ (20) 

The use of (20) for in lends rigor to the thermodynamic definition 
= -d^^da-"'' (Section 1). By extension through V ^ - this 
development also provides a thermodynamic basis for . Thus the center of 
contraction or expansion model makes the notion of formation and relaxation 
volume tensors precise, and also provides a thermodynamic basis for these quan- 
tities. 

Remark 1: Nowick and Heller (1963) have defined a volume tensor analogous 
to . Their expression is vA^j, where v is the atomic volume and Xij is what 
they call an "elastic dipole tensor" . The volume tensor vXij is conjugate to 
the external stress via the interaction energy. However, Xij is "equal to the 
average strain per mole fraction of defects all aligned in a particular orientation" . 
This definition is therefore considerably different from our treatment. To our 
knowledge neither this definition of a volume tensor nor that of Aziz et al. (1991) 
has been rigorously derived in a kinematic or thermodynamic sense. Our search 
of the literature on point defects has not uncovered either a kinematic definition 
of V as in Section 3.1 or a thermodynamic basis for it as in Section 3.2. 

3.3 Atomistic calculations for the formation volume and 
dipole 

The solution for any elastic field around the center of contraction or expan- 
sion depends on the dipole strength, D, which must be known to fully specify 
the elasticity problem. Thus, the formation volume, , cannot be obtained 
from elasticity alone. The dipole, D, and hence , must be obtained from 
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atomistic calculations using empirical interatomic potentials, tight binding or 
ab initio methods. A comparison between such atomistic calculations and a 
simple continuum result allows a determination of D in closed-form as shown 
in Section 3.3.1. With the value of D thus known, continuum elasticity can be 
used for a variety of calculations involving the elastic fields around the point 
defect. In Section 3.4 we show that the elasticity results of Section 3.1 also pro- 
vide important information regarding the validity of atomistic calculations for 
tr[V*]. The results of Section 4.3 similarly demonstrate the importance of in- 
teractions between the vacancy, elastic fields and boundary loading mechanism 
in determining the strain energy of formation, 

We have undertaken atomistic calculations of vacancy formation using the 
Stillinger- Weber interatomic potential (Stillinger and Weber, 1985). Although 
this empirical potential is a rough approximation of silicon, particularly at the 
defect, it is suitable for the purpose of checking consistency between continuum 
elastic and atomistic predictions for and 'W^ in the far-field. 

The Stillinger- Weber potential consists of two- and three-body terms gov- 
erning the interactions of N atoms: 

N N 

$(l,...,iV) = Y.V2ir^i)+ «3(T-„r,,rfc) (21) 

i<j i<j<k 

V2{n,) = (22) 

where energy and length units e and ^ are chosen to give /2 a minimum value 
of —1 if its argument is 2^/^ . The two-body function, /2, depends only on the 
distance = jr^ — rj\ between a pair of atoms with position vectors and 
rj and has a cut-off at ?■ = a without discontinuities in any derivatives with 
respect to r: 

The three-body function, /a, depends on the scalar distances between the three 
atoms and also on the angle subtended at the vertices of the triangle formed by 
the three atoms. In the following relations Ojik is the angle subtended at vertex 
i between atoms with position vectors Tj and rfc: 



f j \ I ^ij Tik \ I Tji Tjk 



h('-f/-^,9.k,] (25) 



Aexp(^:^ + {cos 9 + l)^ -ri and r2 < 



: ri or r2 > a. 
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The equilibrium bond angle in the tetrahedral structure of silicon satisfies 
cos 9 = —1/3. The energetic contribution of any bond angle distortions are thus 
represented by the trigonometric term in (26). The lattice spacing and bond 
energy of silicon at K are obtained with ^ = 0.20951 nm and e = 2.346 eV. 
The optimized set of parameters for this function is: 



Energy minimization calculations were performed using the conjugate gra- 
dient method for system sizes ranging from 64 to 64,000 atoms. The smallest of 
these systems is comparable to moderate-sized ab initio calculations. (We are 
unaware of point defect calculations with ab initio methods having more than 
512 atoms.) The computational cells were cubic in shape. Separate calculations 
were run with with periodic and free boundary conditions respectively. In each 
case the perfect crystal was equilibrated from randomized initial conditions for 
the atoms. Subsequently, the central atom was removed to model a vacancy 
and the system was allowed to equilibrate again. The simulations in periodic 
cells were run at zero average normal traction over each face by including a 
Lagrange multiplier that is energy-conjugate to the cell size. At each iteration 
of the algorithm the total energy was minimized and the cell size was varied to 
obtain zero average normal traction on each of the six faces. 

The energy minimization calculations yielded a mean formation volume over 
the different system sizes of tr[yf] = -13.8 A3 in the minimum energy config- 
uration. This is to be compared with the atomic volume of silicon which is 
f2 = 20 A'^, implying a relaxation volume of tr[V'^] = —33.8 A'^ according to 



Remark 2: The value of tr[y^] = -13.8 A^ suggests a rather strong relaxation 
in silicon, but is typical for the Stillinger- Weber model, which has previously 
been reported to result in large relaxations of the vacancy's nearest-neighbor 
atoms (Balamane et al., 1992). While the quantitative result is not expected to 
be physically accurate, the methodology established in this paper is of greater 
importance. 

3.3.1 Tensorial form of V 

The atomistic calculations also yield the tensorial form of V . The definition 



A = 7.049556277, B 
p 4, g = 0, a 
A = 21.0, 7 



1.80 



1.20. 



0.6022245584 



(16). 




crys 



can be rewritten as 
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where it now is the displacement field of relaxation around the defect obtained 
from the atomistic calculation. Gauss' Theorem then gives 

V-j = j \ {u,n,+u,n{)AA. (27) 

With either periodic or traction-frcc boundary conditions, since the normal, 
n, is fixed for a given surface of the computational cell, it follows that the 
integral in (27) reduces to the symmetric dyadic product of the arcal average 
of each displacement component, J uAA^ and n. The symmetric distribution 
of the components Uj on a face with non-zero normal component n^, for j ^ i 
ensures that V is a diagonal tensor. Symmetry with respect to the three basis 
vectors, {ei, 62, 63}, then ensures isotropy of . These numerical results were 
borne out in the atomistic calculations to a reasonable degree of accuracy. The 
tensors V, obtained for 512 and 64,000 atom calculations using this formula, 
are recorded in Appendix A. The averages and standard deviations over a 
number of samples are reported. 

Equation (27) can be applied to any atomistic calculation to determine Vlj. 
In combination with (16) it uniquely determines the defect dipole tensor, Dij. If 
the atomistic calculation does not involve traction-free boundary conditions, the 
discrete analogue of the surface integral in (16) can be determined very easily 
from the atomistic results by an obvious generalization of the procedure used 
below in Section 3.4. This establishes the thermodynamically-correct relation 
for or V^j. The defect dipole tensor can then be calculated in closed form 
using (16) and (27). 

3.4 Boundary conditions and defect symmetry in deter- 
mining tr[V^] 

In this section wc evaluate the influence of the integral in (16) on the value 
obtained for the scalar relaxation volume, tr[V^]. We consider the cubic com- 
putational cells used in our atomistic calculations. The cell faces were aligned 
with the cubic crystal directions. Wc emphasize that the results of this subsec- 
tion hold where the theory of linear elasticity remains valid. One may therefore 
expect that they will hold in the far-field of the point defect. 

The contribution of the surface integral in (16) needs to be evaluated on 
any one face only since the contributions from the remaining faces are equal by 
symmetry. Consider the [10 0] face, denoted by Si with x = (/, 0:2, xa)"^, where 
/ is the half-length of the computational cell and — / < X2,a;3 < I. The unit 
outward normal to this face is rii = (1,0,0)"^. The components of the traction 
vector = cr^ni, for zero average normal stress, satisfy the following relations 
for an isotropic defect: 

Jt^dA^O, = 0, = Oon 5"!, (28) 



13 



The integral condition on follows from the zero average normal stress con- 
dition. (Hereafter, the phrase "periodic boundary conditions" will imply the 
zero average normal traction condition also, unless explicitly stated otherwise.) 
The pointwise conditions on the shear components ^3 follow since the 

shear stress components (721 and (731 are zero on Si , which in turn follows from 
symmetry across the periodic boundary and isotropy of the defect. Therefore, 
the cubic anisotropy of silicon makes the trace of the boundary integral in (16) 
vanish when evaluated on Si: 



1 



Ciiii + 2C1122 



{t^l + tlx2 + tlxi) dA = (29) 



-I X3- 



by Equation (28). Therefore, according to linear elasticity, the periodic bound- 
ary conditions result in the same scalar relaxation volume as the traction-free 
boundary, provided the defect is isotropic and located at the center of the peri- 
odic cell. From (16) we have 

^'^^'^ = - r TW^ *^[-^] + 

U^llll + /1U1122 

for silicon with cubic anisotropy and an isotropic point defect. Therefore, ac- 
cording to linear elasticity, the calculation of tr[V^] is also unaffected by the 
choice of periodic or traction-free boundary conditions. 

With tr[y] = -33.8 from the atomistic calculation, Cnn = 1.616 x 10" 
Pa and C1122 ~ 0.816 x 10^^ Pa for Stillinger- Weber silicon, (Balamane et al., 
1992) this gives a dipole strength tr[D] = 36.594 x 10~^^ N-m for the isotropic 
vacancy in the Stillinger- Weber model. 

Figure 3 presents atomistic data for tr[V^] = tr[V] + 51 as a function of 
number of atoms used in the calculations. Several calculations were performed 
for each choice of system size (number of atoms). These simulations were started 
from randomized initial positions of atoms. Each data point represents the for- 
mation volume of the calculation with the lowest energy for the given number of 
atoms. The relative independence of tr[V^] obtained by using periodic boundary 
conditions bears out the result in (29). The results from atomistic simulations 
with traction-free boundaries are also shown. The relatively slow convergence of 
these results with system size arises from the fact that the atoms on the surface 
are not identical in bonding coordination to the bulk atoms. This also leads to 
the departure from the elasticity result that tr[V^] is independent of the nature 
of the boundary condition (periodic/free). This discrepancy becomes negligible 
for systems having 64,000 atoms.''' The implication is that independence from 
type of boundary condition is obtained only in the larger systems where elastic 
effects dominate. 



^However, for systems much bigger than 64,000 atoms, the sheer number of atoms intro- 
duces many more spurious minima in the energy minimization calculations. The conjugate 
gradient algorithm tends to get "trapped" in these spurious minima, leading to larger errors 
in the relaxation volumes that are obtained. 
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Figure 3: Comparison of formation volume calculated by atomic relaxation 
for traction-free boundary conditions (Atomistic FB) and periodic boundary 
conditions (Atomistic PB). 

Remark 3: The result of an atomistic calculation for tr[V"^] has been shown to 
be independent of the shape and location of the boundary relative to isotropic 
point defects centered in a computational cell with periodic boundary condi- 
tions. It is therefore independent of the size of the crystal and depends only on 
D, the elasticity of the crystal, and the atomic volume. This type of bound- 
ary condition is therefore well-suited for volume change calculations of isotropic 
defects. 



4 The strain energy of formation 

The observed dependence of on boundary conditions in (16) motivates a 
similar examination of the dependence of the strain energy of formation. In 
the continuum setting it proves convenient to begin such an analysis with an 
infinite crystal and subsequently subject it to cutting and welding operations in 
the manner of Eshelby (1957, 1961) in order to address finite crystals. Let i?crys 
now denote a simply-connected finite subset of the infinite crystal enclosing the 
center of contraction or expansion. Let n denote the unit outward normal to 
i?crys- Let the strain energy in the region i^crys be '^l^^,, and the strain energy 
of the infinite crystal be These energies are related as follows: 



15 



Applying the Divergence Theorem to the integral, and noting that CTj^^- = 
exterior to Bays, leads to 

'^i. = <. - / la°fn,urdV. (30) 

We will not pursue an explicit expression for or ^^^^ . Instead we will 
relate them to the strain energy of a finite crystal with boundary S'crys- 

4.1 Finite crystal with traction- free boundaries 

Now consider turning Bcrys into a finite crystal with traction-free boundaries. 
This can be achieved by holding the traction fixed at cr^rij on S'crys, cutting 
along this surface, and then applying a traction (jfjUj — —a°°nj to S'crys^ where 
afj is the image stress field so produced in i?crys- Let the corresponding dis- 
placement field be denoted by uf , and its strain field by efj = ^ijki'^ki- The 
relation between the strain energy in this configuration, , and '^^^ is 




where we have used the fact that the stress field —crfj applied to the traction-free 
crystal results in a strain field —sjj in ^crysi returning it to the configuration 
with energy "^^^^ . Furthermore, we have used the result that there is no in- 
teraction energy between an external traction, —af^rij, and the internal stress, 
afj + a°°, since the internal stress state is traction free [{afj + cr°j)nj = by 
construction] at the boundary, S'crys- 

Using integration by parts, the Divergence Theorem, and the fact that afj j = 
this is reduced to 

= ^F. + / laf^n.uJdA. (31) 

From (30) and (31) combined with the result (jfjUj = —<T°°nj at S'crys, we 
have 




S'crys 



4.2 Finite crystal of cubic shape and periodic boundaries 

Consider a cube-shaped domain, -Bcrys- Periodic boundary conditions imposed 
on the displacement at constant cell volume imply that the total displacement 
vanishes on Scrys- This state is obtained by subjecting the final configuration 
with traction-free boundaries to the additional boundary displacement uf = 
—u°° — uf at S'crys- Let ef, be the corresponding strain field over Bcrys, and 
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afj = Cijkis]^i be the stress field in Bcrys- Now, the strain energy of the crystal 
is 

where the absence of interaction energy between afj and the internal stress 
with {(tJj + (J°j)'nj = on S'crys has been used. Using integration by parts, the 
Divergence Theorem, and the fact that af^ j — we have 

'^pV =^fV+ / l4^ju[dA. (33) 

5crys 

The implications of (30-33) for the thermodynamics of point defect forma- 
tion are elaborated in Sections 4.3 and 5 below. 

Remark 4: We have used the symbol for the strain energies of formation, 
since, for the purely mechanical processes considered in this paper, the change 
in internal energy, , is the strain energy when we only consider the crystal 
away from the defect. An additional contribution to the formation energy arises 
from the change in bonding at the defect, which is, of course, not represented 
in elasticity. 

4.3 Numerical evaluation of strain energies of formation 

The strain energy of formation with traction-free boundaries on a cubic com- 
putational cell was numerically-evaluated as follows: The stress at the bound- 
aries of the cubic computational cell embedded in an infinite crystal was ob- 
tained from (10) and a°° = Cijkis'^. The resulting force distribution on a 
surface, Sa , of the computational cell was considered, where the surface normal 
was denoted by = j^e|Q|, a ~ {±1,±2,±3}. The subscript {•)\a\ de- 
notes (•)!, (•)2 or (•)3. The force distribution was represented by point forces 
^jm(^a/) = cF'ijixM)najAM, with Am, M = 1, . . ■ N^^ being the area associ- 
ated with N^^ points, xm, on each surface of the computational cell. These 
points subsequently defined the surface nodes of a finite element discretization 
of the computational cell with N^^^ nodes. The finite element mesh had equal- 
sized, trilinear, cubic elements. For such a mesh. Am equals the area of each 
face of the cubic elements for a point that docs not lie on an edge or vertex 
of the computational cell; Am equals half the area of each face of the cubic 
elements for a point on an edge; and Am equals quarter the area of each face of 
the cubic elements for a point on a vertex. 

Image forces were applied on the nodes of Sa by reversing the above-computed 
nodal forces so that, when superposed on the tractions a![jna-, the resulting 
field at the boundary nodes was traction-free. The displacement field uf due 
to these nodal image forces on the boundary was obtained as the finite element 
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displacement solution, and the integral in (32) was calculated as 



. 3 N,i, 



a— —3 1 



(34) 

This is the difference in energies between the infinite crystal and the traction- 
free boundary case. The calculations were carried out for tr[D] = 36.594 x 10~^^ 
N-m. 

The difference between strain energies of the periodic boundary and free 
boundary cases, , — given by the integral in (33). Recall that in the 

atomistic calculations with periodic boundary conditions, the computational cell 
was also relaxed to have zero average normal traction, and pointwise zero shear 
tractions on the surfaces. In the finite element calculations this was achieved 
by applying the following traction field over the surface Sa'. 



— - ' (35) 



Arca(5„) / (-4"". ) ' * = 



a 



-af.na- , Hi ^ \a\ 



Note that while tf is constant over the surface Sa for i = jal, it varies over Sa 
for i \a\. 

Let the displacement field result from application of traction tf' defined 
as in (35) over each surface, Sa- Since tf is constant over Sa for i = \a\, and 
the shear components do not give rise to normal strains in a cubic material, it 
follows that uf is constant over Sa for i = \a\. 

In order to extend (33) to the periodic boundary condition calculation with 
zero average normal traction, the integral in that equation was replaced by 

3 



a— — 3c 



where — ~ — uf . From (35) and uf = constt on Sa for i = \a\, this 
integral is 



a=-3o 

The finite element displacement solution obtained for —uf was added to the 
infinite crystal displacement field —u°°- This field,— — wf , evaluated on the 
boundaries, Sa, was then reapplied as a displacement boundary condition to 
the finite element mesh to obtain the surface traction field afjUa^ Then (35) 
and (36) gave the required energy difference. 

Figure 4 is a comparison of the energy differences — '^/^ and — '^^^ 
obtained with the atomistic and elasticity calculations. The actual value of 
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has not been obtained directly from elasticity. However, the integrals in 
(34) and (36) both are 0(r~^), where 2r is the characteristic linear dimension 
of the computational cell. Therefore, and '^p^ converge to as r — > 

oo. Since cx A^~^, where is the number of atoms in the system, this 
convergence is confirmed in Figure 4. Analysis of (32) and (33), or of (34) 
and (36) indicates that '^p^ and '^^v converge to from below. Figure 4 
confirms this prediction, since the energy differences — '^l^,-^ are all positive. 
The trends are replicated by the atomistic calculations, in which 'i^-^^. and '^/^ 
are directly calculated. The asymptotic value to which and '^p^ converge 
in the atomistic calculations has been used as . 

The discrepancy between linear elasticity and atomistic calculations for the 
smaller cells (e.g. 512 atoms) in Figure 4 is most probably due to the highly 
nonlinear deformations around the vacancy that dominate at these small scales. 
The 512-atom computational cell is made up of 4 x 4 x 4 Si unit cells. The 
boundary is only 10.86 A away from the vacancy. At such small scales linear 
elasticity provides a poor estimate of the strain energy of the highly distorted 
lattice around the vacancy. 

The difference between the energies of the traction-free and periodic bound- 
aries with the atomistic calculations in Figure 4 most probably is a manifestation 
of an interaction energy in the near field of the defect: The periodic boundary 
configuration effectively applies a boundary traction to the traction-free config- 
uration, as explained above. This external stress field interacts with the internal 
stress of the traction-free boundary configuration. Such an interaction energy, 
however, vanishes according to the theory of linear elasticity. This result has 
been employed in Sections 4.1 and 4.2, and is possibly responsible for the very 
small differences observed between '^p^ and '^p^ using linear elasticity. For the 
smaller cells, the strong nonlinear effect possibly invalidates this result of lin- 
ear elasticity causing the discrepancy in '^p^ — '^/^ between the atomistic and 
linear elasticity calculations. The trends, however, are the same: The free and 
periodic boundary cases converge from below to the infinite crystal result. The 
periodic boundary case is higher in strain energy due to the imposed constraint 
that preserves the flat shape of the cubic cell's boundaries. This difference, 
while small, exists in the elasticity calculations also as seen in Table 1. 

5 Conclusions 

The central conclusions from this work are the following: 

1. Equation (16) implies that the formation volume of a point defect is a 
well-defined quantity. It is to be measured on a crystal whose surfaces 
are traction-frcc. Defined in this manner, it depends only on the strength 
of the defect, the elasticity of the material and the atomic volume. It 
is independent of the crystal's shape and size, and of the location of the 
defect within the crystal. 
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Table 1: Comparison of strain energy of vacancy formation in a cubic cell with 
periodic and traction-free boundaries — elasticity results. For comparison, the 
atomistic calculations converge to = 2.8237 eV. Note that the cell sizes, 
expressed as number of Si atoms that make up a crystal of the same size, do not 
correspond exactly to the cell sizes actually used in the atomistic calculations. 



Cell size in number of atoms 


^py - ^Fv (^eV) 


512 


45.791 


4096 


5.759 


13824 


1.714 


46656 


0.510 


110592 


0.215 



2. On the basis of Item 1 above. Equation (16) provides an exact quantitative 
measure of the appropriateness of any class of boundary conditions used 
on an atomistic calculation to extract formation volumes of point defects. 

3. The derivation ending in Equations (19) and (20) demonstrates that Es- 
helby's calculation (17) on the work done by external stress on a trans- 
forming inclusion has an analogue for point defects. This is an exact result. 
It provides a rigorous thermodynamic basis for the concepts of formation 
and relaxation volume tensors that goes beyond formal definitions of the 
type = ~a^V9cr'=>=. 

4. The combination of (27) and (16) provides a closed-form expression for the 
defect dipole tensor that can be determined from any atomistic calculation. 

5. Equations (30-33) quantify the influence of traction-frcc and periodic 
boundary conditions on the strain energy of formation of point defects 
according to continuum linear elasticity. In each case, the corresponding 
boundary integrals can be calculated. Significant discrepancies are ob- 
served between the continuum clastic and atomistic results for the smallest 
system sizes (512 atoms) because the highly nonlinear lattice distortion 
around the defect permeates the entire computational cell when the lat- 
ter is composed of a small number of atoms. Linear elasticity is not an 
accurate theory for such nonlinear distortions. Additionally, for atomistic 
systems of this size, there appears to be an interaction energy between 
an externally-applied boundary traction and an internal stress field with 
vanishing traction at the boundary. According to linear elasticity this in- 
teraction energy vanishes exactly. This could lead to the much smaller 
diff'erence in strain energies of formation between the traction-free and 
periodic boundary cases obtained via linear elasticity. Indeed, this sug- 
gests that (31-33) are not accurate for systems smaller than roughly 4000 
atoms. 

6. The difference between strain energies of vacancy formation obtained with 
traction-free boundaries and periodic boundaries (that are traction free 
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Figure 4: Comparison of strain energies of formation calculated by atomic re- 
laxation and elasticity for periodic (PB) and traction-free (FB) boundary con- 
ditions. 



on average) vanishes with increasing cell size. Linear elasticity predicts 
that the scalar vacancy formation volumes, however, are identical for an 
isotropic defect that is centered in the computational cell. Yet, the atom- 
istic results show significant differences for the scalar formation volume 
for the two types of boundary conditions (Figure 3) at smaller cells. This 
discrepancy can be attributed to the larger influence of boundary atoms' 
coordination numbers for small cells, and it diminishes for the larger calcu- 
lations. We expect that this discrepancy will be enhanced, due to elastic- 
ity effects, for the much larger class of anisotropic point defects: vacancies 
with the Jahn- Teller distortion and most interstitial configurations, as well 
as for defects not located at the center of the periodic cell. 



A Relaxation volume tensors from atomistic cal- 
culations 

Representative relaxation volume tensors obtained from the energy minimiza- 
tion calculations on systems with varying numbers of atoms, and with traction- 
free or periodic boundary conditions are reported below. For the sake of brevity 
values are reported for only the 512- and 64,000-atom systems. Averages and 
standard deviations (SD) were calculated by running several calculations with 
starting positions of atoms varying randomly from the equilibrium positions in a 
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perfect crystal. The larger standard deviations for larger systems are associated 
with the spurious minima discussed in Footnote 4. 
Traction-free boundary conditions with 512 atoms: 



-13.676 
0.683 
0.695 



0.683 
-13.669 
0.678 



0.695 
0.678 

--13.674 



SD = 



0.192 
0.030 
0.031 



0.030 0.031 
0.172 0.027 
0.027 0.151 



Periodic (traction-frcc on average) boundary conditions with 512 atoms: 



-12.009 
0.706 
0.706 



0.706 
-12.009 
0.706 



0.706 
0.706 
-12.009 



SD = 



0.001 
0.002 
0.002 



0.002 0.002 
0.001 0.001 
0.001 0.001 



Traction-frcc boundary conditions with 64,000 atoms: 



V = 



-11.625 
-0.022 
-0.021 



-0.022 

-11.484 

-0.029 



-0.021 
-0.029 
-11.421 



A^ 



SD = 



0.535 
0.178 
0.185 



0.178 0.185 
0.308 0.138 
0.138 0.339 



Periodic (traction-frcc on average) boundary conditions with 64,000 atoms: 



-11.571 
-0.147 
-0.130 



-0.147 
-11.217 
-0.054 



-0.130 
-0.054 
-11.136 



A^ 



SD 



0.653 0.115 0.057 
0.115 0.473 0.142 
0.057 0.142 0.363 
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